Occupancy Assessment for Bats of the Pacific Northwest

Preliminary Results

Authors
Affiliations

Trent VanHawkins

On Belay Statistical Consulting

Beth Ward

Oregon State University-Cascades

Tom Rodhouse, PhD

National Park Service

Code
library(tidyverse)
library(here)
library(latex2exp)
library(sf)
library(rjags)
library(coda)
library(gt)
library(kableExtra)

setwd("/Users/trentvanhawkins/Library/Mobile Documents/com~apple~CloudDocs/Documents/OnBelay/Projects/BatHub/BatHub_Trend")

# Read in Data
## Full 
dir <- "DataProcessed/results/jags/full/fits/"

filenames <- list.files(here(dir))

#read in the files into a list object
occ_jags <- lapply(filenames,function (x) readRDS(here(paste0(dir, x))))

# rename the list elements to have the correct names
names(occ_jags) <- str_split_i(filenames, "_", 1)

## Sensitivity
dir.sens <- "DataProcessed/results/jags/ORWA_Only/fits/"

filenames.sens <- list.files(here(dir.sens))

#read in the files into a list object
occ_jags_sens <- lapply(filenames.sens,function (x) readRDS(here(paste0(dir.sens, x))))

# rename the list elements to have the correct names
names(occ_jags_sens) <- str_split_i(filenames.sens, "_", 1)

## define some species we will exclude
exclude <- c("tabr")

Introduction

This purpose of this report is to document an occupancy trend analysis for bat species of the Pacific Northwest. Previously, the Northwestern Bat Hub for Population Monitoring and Research (NW Bat Hub) and its partners have documented evidence of a region-wide decline for Hoary Bats between 2016 and 2018 (T. J. Rodhouse et al. 2019). Since publication of this research, the bat hub has continued population research and monitoring under North American Bat Program (NABat) summer-monitoring protocols and has since grown to include surveys in all of Oregon, Washington, and Idaho. Here, we present the preliminary results of an updated analysis including data from 2016-2022 for all of the NABat sample units encompassed by Oregon, Washington, and Idaho.

Methods

Methods for the current study are similar to those published in Rodhouse et. al, 2019 and Kéry and Royle (2021). We will provide brief detail on important components of the data and modeling procedures but further details can be found in the literature.

The Data

We queried data from NW Bat Hub database files for all included up to the time of this report (2016-2022). The NABat summer sampling design and its explanation are well-documented in the literature (Reichert et al. 2021; T. Rodhouse and Irvine 2017; T. J. Rodhouse et al. 2019; Udell et al., n.d.).

Briefly, NABat uses a spatial grid covering the continental United States in 10x10 km sample units—hereto referred to as “grid cells” or simply “cells”— with an ordered sampling priority to ensure a spatially balanced sample that meets statistical modeling assumptions. The current monitoring protocol utilized by the NW Bat Hub calls for one night of acoustic monitoring at each of four spatial replicates within a cell during the summer season. Ideally, each of the spatial replicates are distributed to one of the inter-cardinal 5x5 km quadrants (NE, NW, SE, SW) of a cell, but when this is not possible multiple detectors may be deployed in a single quadrant. The current modeling strategy is also flexible to >4 spatial replicates, but is not required. In instances where detectors were deployed for >1 night, only the first night was subsetted for analysis.

The Mexican Free-tailed Bat (Tadarida brasiliensis) has been removed from this analysis due to small sample size and limited range.

The Model

The Bayesian framework used for this analysis is thoroughly described in T. J. Rodhouse et al. (2019) and Kéry and Royle (2021), but we will briefly define the model and its parameters here. The framework is broken down into two sub-models, an occupancy sub-model and a detection sub-model.

First, we have a true underlying occupancy status \(Z_{it}\) at site \(i\) and year \(t\). This random variable can take on the value of 1 (occupied) or 0 (unoccupied). The true occupancy status is unknowable (since we cannot truly know that a site is unoccupied via acoustic monitoring), but we can say that it arises from a probability distribution with some probability that the site is occupied (\(\psi_{it}\)).

\[Z_{it} \sim \text{Bernoulli}(\psi_{it}),\]

Furthermore, we assume that the parameter \(\psi_{it}\) can be partially explained by a set of environmental covariates, \(\bf{x}\) (Table 2), for which we will estimate an effect \(\alpha\). We will also allow the probability of occupancy to vary between years by explicitly modeling the probabilities of transitioning occupancy states in years \(\{2, ..., T\}\). The probability of colonization is denoted by \(\gamma_{t}\) and describes the probability that a an unoccupied (\(Z_{i(t-1)} = 0\)) cell in year \(t-1\) transitions to being occupied (\(Z_{i(t)} = 1\)) in year \(t\). Conversely, The probability of survival is denoted by \(\phi_{t}\) and describes the probability that a an occupied (\(Z_{i(t-1)} = 0\)) cell in year \(t-1\) remains occupied (\(Z_{i(t)} = 1\)) in year \(t\). This model parameterization is termed “fully dynamic” by Kéry and Royle (2021).
\[\text{logit}(\psi_{i1}) = \alpha_{01} + {\bf x}_i\alpha,\] \[\text{logit}(\psi_{it}) = \gamma_{t} + {\bf x}_i\alpha + Z_{i(t-1)}\phi_{t}\text{ for }t \in \{2,\ldots,T\},\]

Finally, we have a detection sub-model which estimates the conditional probability that we detected a bat species (\(p_{ijt}\)) at site \(i\) and replicate \(j\) on year \(t\) given that site \(i\) was occupied at year \(t\). We assume that \(p_{ijt}\) arises from a similar probability distribution and may also be explained by some environmental covariates \(\bf{v}\), for which we wish to estimate effects (\(\bf{\beta}\)).
\[[Y_{ijt} \mid Z_{it}=1] \sim \text{Bernoulli}(p_{ijt})\] \[\text{logit}(p_{ijt}) = \beta_0 + {\bf v}_{ijt}\beta,\]

A full summary of these parameters is provided in Table 1 and environmental covariates in Table 2. For the purpose of this report, we are most interested in reviewing the posterior probability of occupancy (\(\psi_{it}\)), effect estimates for occupancy environmental covariates (\(\bf{\alpha}\)), and effect estimates for detection environmental covariates \(\bf{\beta}\). MCMC sampling was conducted in 4 chains, each run for 5000 iterations with a thinning rate of 3 and burn-in of 1000 iterations. Vague priors (Normal (mean = 0, var = 10)) priors were used universally. All models were fit using rjags ver. 4.16. Code used to produce this analysis can be found in Section 6.

Estimating Trend

The goal of this analysis is to estimate an occupancy trend for each of the species included in this analysis. There are many possible ways to accomplish this goal, but the one used in T. J. Rodhouse et al. (2019) is to define a variable

\[ \hat\lambda = \frac{\hat\psi_{(2022)}}{\hat\psi_{(2016)}} \]

This comparison of occupancy probabilities allows estimation of relative occupancy over the duration of this study, but is no substitute for examining the posterior distributions of \(\hat\psi_{it}\).

Table 1: Model parameters and descriptions.
Parameter Description
\(\psi_{it}\) Posterior occupancy probability for site \(i\) at time \(t\)
\(\alpha_{01}\) Occupancy model intercept at year \(t = 1\)
\(\gamma_{t}\) Probability of colonization at year \(t = 2, \ldots, T\)
\(\bf{\alpha}\) Occupancy model effect size for environmental covariates
\(\bf{x}\) Occupancy model environmental covariates Table 2 (a).
\(\phi_{t}\) Autoregressive parameter allowing for different occupancy probabilities on year \(t = 2,\ldots, T\)
\(p_{ijt}\) Conditional probability of species detection at site \(i\) and replicate \(j\) on year \(t\) given site \(i\) is occupied at year \(t\)
\(\beta_0\) Detection model intercept
\(\bf{v}\) Detection model environmental covariates Table 2 (b)
\(\bf{\beta}\) Detection model effect sizes for environmental covariates
Table 2: Covariates included in (A) the occupancy model, or (B) the detection model. All covariates were first standardized to mean = 0 and sd = 1 to promote computational efficiency with the exception of Clutter (centered only) and categorical Waterbody. Forest Cover and Cliff Cover were natural-log transformed before standardization.
(a) Environmental covariates included in the occupancy sub-model.
Name Description
Forest Cover (%) Percentage of the sample unit that classified as forest cover (NLCD)
Precipitation (mm) Mean annual precipitation (WorldClim)
Elevation (m) Maximum Elevation (DEM)
Cliff Cover (%)

Percentage of the sample unit classified as cliff cover (LandFire)

(Anpa, Euma, Myci, Pahe)

(b) Environmental covariates included in the detection sub-model.
Name Description
Minimum Temperature Minimum temperature (\(C^{\circ}\)) on the night of detector deployment
Day Length Day length (previous day) on the night of deployment
Clutter (%) Categorical clutter (0 = 0%; 1 = 1-25%; 2 = 26-50%; 3 = 51-75%; 4 = 76-100%)
Waterbody Presence of a waterbody at deployment site (1 = yes; 0 = no)

Sensitivity Analysis

Previous studies released by the NW Bat Hub and its partners have focused on acoustic records captured in the 2016 to 2019 field seasons and, therefore, exclude Idaho. In order to compare model results from this analysis with those more comparable to previous research, we have conducted a sensitivity analysis using data from only Oregon and Washington. Methods for this sensitivity analysis are the same as for the primary analysis and will be presented alongside results from the full analysis.

Results

Code
occ_data <- readRDS(here("DataProcessed/results/jags/occ_data.rds"))
occ_data_sens <- readRDS(here("DataProcessed/results/jags/occ_data_sens.rds"))

n_sites <- occ_data$laci$n_sites
n_sites_sens <- occ_data_sens$laci$n_sites

A total of 355 independent cells were included across the seven year period from 2016 to 2022 in the primary analysis and 254 in the sensitivity analysis (Table 3; Figure 1).

Code
occ_map <- readRDS(here("DataProcessed/results/stan/full/occ_map.rds"))

occ_map$anpa$means_map_long %>% 
  st_drop_geometry() %>% 
  filter(sampled == 1) %>% 
  group_by(years) %>% 
  reframe("Sample Units" = sum(sampled)) %>% 
  rename("Year" = years) %>% 
  gt() %>% 
  tab_options(table.width = pct(50)) %>% 
  cols_align("left")
Table 3: Sample units surveyed by year
Year Sample Units
2016 82
2017 60
2018 145
2019 174
2020 242
2021 267
2022 294

Partners at in Idaho joined the NABat effort coordinated by the NW Bat Hub during the 2020 season, and sampling efforts are reflected in Figure 1.

Code
occ_map$anpa$means_map_long %>%
  mutate(sampled = factor(sampled, labels = c("No", "Yes"))) %>% 
  ggplot()+
  geom_sf(aes(fill = factor(sampled)))+
  facet_wrap(~years)+
  scale_fill_discrete(type = c("transparent", "blue"))+
  labs(fill = "Sampled")+
  theme_bw()+
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.x = element_blank(),
          axis.ticks.y = element_blank(),
        legend.position = "bottom")
Figure 1: Map of NABat grid cells (10x10 km) sampled by year.

Modeling Results

All models satisfied MCMC convergence criteria with nominal \(\hat{R}\) (<1.1) and Effective Sample Size (>500) values. Adequate chain mixing was assessed by visual inspection of traceplots.

Covariates

Figure 2 gives summaries of the posterior distributions for each environmental covariate from the occupancy sub-model, for each species included in the analysis.

Code
# Extract mcmc draws and convert to dataframe for summarization and plotting (outputs all.jags which contains all mcmc draws of all params for all species)
mcmc_extract <- function(occ_jags){
  for(i in 1:length(occ_jags)){
    cur.fit <- occ_jags[[i]]
    spp <- names(occ_jags)[i]
  
    for(j in 1:length(cur.fit)){
      tmp <- cur.fit[[j]] %>% 
        as.data.frame() %>% 
        mutate(iter = seq(1:nrow(cur.fit[[j]])),
               chain = j,
               spp = spp)
      if(j == 1){dat <- tmp}
      else{dat <- bind_rows(dat, tmp)}
    }
    if(i == 1){all.jags <- dat}
    else{all.jags <- bind_rows(all.jags, dat)}
  }
  return(all.jags)
}
## Use our function to extract our data
all.jags <- bind_rows(mcmc_extract(occ_jags) %>% mutate(analysis = "OR|WA|ID"), mcmc_extract(occ_jags_sens) %>% mutate(analysis = "OR|WA Only"))


## Summarise parameter estimates
param_summ <- all.jags %>%
  select(dplyr::contains('alpha') | dplyr::contains('beta'), spp, analysis) %>% 
  rename("Min. Temp" = beta1,
         "Daylight" = beta2,
         "Clutter 0" = beta3,
         "Clutter 1" = beta4,
         "Clutter 2" = beta5,
         "Clutter 3" = beta6,
         "Water" = beta7,
         "Log(Forest Cover)" = `alphas[1]`,
         "Precipitation" = `alphas[2]`,
         "Elevation" = `alphas[3]`,
         "Log(Cliff Cover)" = `alphas[4]`) %>% 
  pivot_longer(cols = -c(spp, analysis), names_to = "param") %>% 
  group_by(spp, analysis, param) %>% 
  summarise(mean = mean(value),
            q2.5 = quantile(value, probs = c(0.025), na.rm = T),
            q25 = quantile(value, probs = c(0.25), na.rm = T),
            q50 = quantile(value, probs = c(0.5), na.rm = T),
            q75 = quantile(value, probs = c(0.75), na.rm = T),
            q97.5 = quantile(value, probs = c(0.975), na.rm = T)) %>% 
  ungroup() %>% 
  mutate(type = if_else(param %in% c("Log(Forest Cover)", "Precipitation", "Elevation", "Log(Cliff Cover)", "alpha01"), "Occurrence", "Detection")) 
  


##Plot alphas
param_summ %>% 
  filter(type == "Occurrence") %>% 
ggplot(aes(x = spp, y = mean, group = analysis))+
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, size = 0.75,
                position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = q25, ymax = q75, color = analysis), width = 0,
                size = 1.5, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  facet_grid(param ~ ., scales = 'free') +
  geom_hline(yintercept = 0, lty = 2) +
  xlab('Species') + 
  ylab('Posterior Distribution (Log-Odds)') +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  ggtitle('Comparing occupancy coefficients (Mean & 95% CI)',
          'Multi-season, single-species model')+
  labs(color = "Legend")
Figure 2: Summary of posterior distributions for occupancy model covariates, by species. Point locations represent the posterior means, colored regions represent the 50% credible interval, and black lines represent the 95% credible interval. Outputs provided on the log-odds scale.

Figure 2 shows stable estimates between the full and sensitivity analyses. Estimates are reported on the log-odds scale and can be generally be interpreted as “a change of 1 standard deviation in X is associated with an increase of Y in the log-odds of occupancy”.

Figure 3 gives summaries of the posterior distributions for each environmental covariate from the detection sub-model, for each species included in the analysis.

Code
param_summ %>% 
  filter(type == "Detection") %>% 
ggplot(aes(x = spp, y = mean, group = analysis))+
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, size = 0.75,
                position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = q25, ymax = q75, color = analysis), width = 0,
                size = 1.5, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  facet_grid(param ~ ., scales = 'free') +
  geom_hline(yintercept = 0, lty = 2) +
  xlab('Species') + 
  ylab('Posterior Distribution (Log-Odds)') +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  ggtitle('Comparing detection coefficients (Mean & 95% CI)',
          'Multi-season, single-species model')+
  labs(color = "Legend")
Figure 3: Summary of posterior distributions for detection model covariates, by species. Point locations represent the posterior means, colored regions represent the 50% credible interval, and black lines represent the 95% credible interval. Outputs provided on the log-odds scale.

Figure 3 again shows stable estimates for model coefficients across species. These results may also be interpreted as “a change of 1 standard deviation in X is associated with an increase of Y in the log-odds of detection, given that the species is present.”

Occupancy

Figure 4 gives a summary of the posterior distribution of estimated occupancy probabilities across all predicted grid cells in the three-state region of Oregon, Washington, and Idaho. Again, these are preliminary results and are subject to change. Exact estimates and credible intervals are available in Table 4 and Table 5 for sensitivity results.

Code
psi_summ <- all.jags %>% 
  select(contains("psi"), spp, analysis) %>% 
  pivot_longer(-c(spp, analysis)) %>% 
  mutate(year = case_when(str_detect(name, '1')~2016,
                          str_detect(name, '2')~2017,
                          str_detect(name, '3')~2018,
                          str_detect(name, '4')~2019,
                          str_detect(name, '5')~2020,
                          str_detect(name, '6')~2021,
                          str_detect(name, '7')~2022),
         type = if_else(str_detect(name, 'avg'), 'Average', 'Hat')) %>% 
  group_by(spp, analysis, type, year) %>% 
  summarise(mean = mean(value),
            q2.5 = quantile(value, probs = c(0.025), na.rm = T),
            q25 = quantile(value, probs = c(0.25), na.rm = T),
            q50 = quantile(value, probs = c(0.5), na.rm = T),
            q75 = quantile(value, probs = c(0.75), na.rm = T),
            q97.5 = quantile(value, probs = c(0.975), na.rm = T)) %>% 
  ungroup()

psi_summ %>% 
  filter(type == "Hat") %>% 
  ggplot(aes(x = as.factor(year), y = mean, group = analysis))+
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, size = 0.75,
                position = position_dodge(width = 0.25)) +
  geom_errorbar(aes(ymin = q25, ymax = q75, color = analysis), width = 0,
                size = 1.5, position = position_dodge(width = 0.25)) +
  geom_point(position = position_dodge(width = 0.25)) +
  theme_bw() +
  facet_grid(spp ~ ., scales = 'free') +
  xlab('Year') + 
  ylab(expression('Posterior Distribution' ~ hat(psi)[t])) +
  ylim(c(0,1))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  ggtitle('Occurrence Probabilities Over Time',
          'Multi-season, single-species model')+
  labs(color = "Legend") 
Figure 4: Region-wide summary of the posterior distribution of \(\hat{\psi_t}\) along with 50% credible interval (colored region) and 95% credible interval (black line).

Trend Estimation

As mentioned previously, we define \(\lambda = \frac{\hat\psi_{2022}}{\hat\psi_{2016}}\). This provides a relative comparison of trend over the endpoints of our study period.

Code
lambda_summ <- all.jags %>% 
  select(contains("lam.tot"), spp, analysis) %>% 
  pivot_longer(-c(spp, analysis)) %>% 
  group_by(spp, analysis) %>% 
  summarise(mean = mean(value),
            q2.5 = quantile(value, probs = c(0.025), na.rm = T),
            q25 = quantile(value, probs = c(0.25), na.rm = T),
            q50 = quantile(value, probs = c(0.5), na.rm = T),
            q75 = quantile(value, probs = c(0.75), na.rm = T),
            q97.5 = quantile(value, probs = c(0.975), na.rm = T)) %>% 
  ungroup()

lambda_summ %>% 
  ggplot(aes(x = spp, y = mean, group = analysis))+
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, size = 0.75,
                position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = q25, ymax = q75, color = analysis), width = 0,
                size = 1.5, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  xlab('Year') + 
  ylab(expression('Posterior Distribution' ~ hat(lambda))) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  ggtitle('Relative Occupancy Probability in 2022 vs. 2016')+
  labs(color = "Legend") 
Figure 5: \(\hat{\lambda}\) estimated for each species.

Estimates for coto and euma in Figure 5 appear to be extremely more likely to occur in 2022 than 2016. This is most likely due to the fact that they have very low occurrence probabilities in 2016, so even small gains in improvement would lead to large values of \(\lambda\). We review the figure without those species in Figure 6.

Code
lambda_summ %>% filter(!spp %in% c("coto", "euma")) %>% 
  ggplot(aes(x = spp, y = mean, group = analysis))+
  geom_errorbar(aes(ymin = q2.5, ymax = q97.5), width = 0, size = 0.75,
                position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = q25, ymax = q75, color = analysis), width = 0,
                size = 1.5, position = position_dodge(width = 0.5)) +
  geom_point(position = position_dodge(width = 0.5)) +
  theme_bw() +
  geom_hline(yintercept = 1, lty = "dashed")+
  xlab('Year') + 
  ylab(expression('Posterior Distribution' ~ hat(lambda))) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  ggtitle('Relative Occupancy Probability in 2022 vs. 2016')+
  labs(color = "Legend") 
Figure 6: \(\hat{\lambda}\) estimated for each species with coto and euma removed to better visualize results. The dashed horizontal line indicates \(\hat{\lambda} = 1\). 95% credible intervals falling strictly below this line indicate a species that is less likely to occur in 2022 than 2016.

We see a decline in relative occurrence probabilities (\(\hat\lambda < 1\)) between 2016 and 2022 for EPFU, MYCA, MYEV, MYTH, and MYVO. Evidence is particularly strong for EPFU, MYTH, and MYVO.

Occupancy maps

Maps of the posterior predictions for all species displayed in Figure 4 are included with this report and visualize the mean and 95% credible intervals of the posterior distribution of \(\hat{\psi_t}\). Figure 7 gives an example of the map for mean posterior occupancy probabilities and Figure 8 displays uncertainty around those means.

Code
psi.post.map <- readRDS(here("DataProcessed/results/jags/psi_post_map.rds"))

psi.post.map %>% 
  filter(spp == "myci") %>% 
  ggplot()+
    geom_sf(aes(fill = mean), size = 0.1, lwd = 0.05) +
    facet_wrap(. ~ year, ncol = 4) +
    viridis::scale_fill_viridis(limits = c(0, 1)) +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.x = element_blank(),
          axis.ticks.y = element_blank(),
          legend.position = "bottom")+
    ggtitle('MYCI Posterior Mean Occupancy',
            'Multi-season, single-species model')+
  labs(fill = "Posterior Mean")
Figure 7: Mean posterior occupancy probabilities visualized for the small-footed bat. Cool colors indicate low probability of occurrence on year t and warm colors indicate a high probability.
Code
psi.post.map %>% 
  filter(spp == "myci") %>% 
    ggplot()+
    geom_sf(aes(fill = width), size = 0.1, lwd = 0.05) +
    facet_wrap(. ~ year, ncol = 4) +
    viridis::scale_fill_viridis(limits = c(0, 1), option = "magma") +
    theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.x = element_blank(),
          axis.ticks.y = element_blank(),
          legend.position = "bottom")+
    ggtitle('MYCI Posterior 95% Credible Interval (Width)',
            'Multi-season, single-species model')+
  labs(fill = "95% CI")
Figure 8: 95% credible interval of posterior occupancy probabilities visualized for the small-footed bat. Warm colors indicate higher uncertainty.

Discussion

This report represents a significant update to monitoring of regional trends for bat species of the pacific northwest. We used the fully dynamic occupancy model proposed by Kéry and Royle (2021) and implemented by T. J. Rodhouse et al. (2019) to summarize occupancy probabilities in the three-state region of Oregon, Washington, and Idaho for each of 14 bat species. We note a measurable decline in occupancy probabilities for Eptesicus fuscus, Myotis volans, and Myotis thysanodes across the three state region of Oregon, Washington, and Idaho.

Dynamic occupancy models have proven to be a powerful tool region-wide population monitoring as they can directly model the latent ecological processes of occupancy, colonization, survival, and extinction and provide optimal estimates for occupancy probabilities (\(\hat\psi\)) across years. Future work may want to focus on developing more precise tools for linear and non-linear estimates of the effect of time (in this case year) on the log-odds of occupancy, as these parameters could provide a clean and clear interpretation of these models.

Acknowledgements

The Northwest Bat Hub thanks its contributing partners for their continued support to better understand regional status and trends of bats in the Pacific Northwest.

Appendix

Figures & Tables

Code
library(readxl)
wbwg <- read_xlsx("/Users/trentvanhawkins/Library/Mobile Documents/com~apple~CloudDocs/Documents/OnBelay/Projects/BatHub/BatHub_SDM/Background/SpeciesNatHistMatrix.xlsx", sheet = 2)

name.match <- wbwg %>% select(fourlettername, SPECIES) %>% drop_na()
psi_summ <- psi_summ %>% left_join(name.match, by = c("spp" = 'fourlettername'))
psi_summ %>% 
  filter(type == "Hat", analysis == "OR|WA|ID") %>% 
  select(SPECIES, year, mean, q2.5, q25, q75, q97.5) %>% 
  mutate(year = as.factor(year)) %>% 
  gt(groupname_col = "SPECIES") %>%
  cols_label("mean" ~ "Mean",
             "year" ~ "Year",
             "q2.5" ~ "2.5%",
             "q25" ~ "25%",
             "q75" ~ "75%",
             "q97.5" ~ "97.5%") %>% 
  tab_spanner(label = "Quantile", columns = 4:7) %>% 
  fmt_number(decimals = 2) %>% 
  tab_options(row_group.as_column = T)
Table 4: Region-wide summary of the posterior distribution of \(\hat{\psi_t}\) with credible intervals for the analysis including Oregon, Washington, and Idaho
Year Mean Quantile
2.5% 25% 75% 97.5%
Antrozous pallidus 2016 0.22 0.11 0.17 0.27 0.39
2017 0.23 0.12 0.18 0.26 0.38
2018 0.17 0.08 0.14 0.20 0.26
2019 0.15 0.08 0.12 0.18 0.24
2020 0.12 0.07 0.10 0.14 0.19
2021 0.14 0.08 0.12 0.16 0.21
2022 0.14 0.08 0.12 0.16 0.21
Corynorhinus townsendii 2016 0.08 0.02 0.05 0.11 0.19
2017 0.14 0.05 0.10 0.17 0.27
2018 0.15 0.07 0.12 0.18 0.26
2019 0.17 0.09 0.14 0.20 0.27
2020 0.18 0.11 0.15 0.20 0.26
2021 0.18 0.12 0.15 0.20 0.25
2022 0.15 0.09 0.12 0.17 0.22
Eptesicus fuscus 2016 0.82 0.70 0.78 0.85 0.91
2017 0.83 0.73 0.80 0.86 0.91
2018 0.82 0.74 0.79 0.84 0.89
2019 0.67 0.57 0.64 0.70 0.76
2020 0.72 0.65 0.70 0.75 0.79
2021 0.69 0.62 0.67 0.72 0.77
2022 0.64 0.57 0.62 0.67 0.72
Euderma maculatum 2016 0.03 0.00 0.01 0.04 0.07
2017 0.08 0.04 0.06 0.09 0.14
2018 0.07 0.04 0.06 0.08 0.12
2019 0.07 0.03 0.05 0.08 0.10
2020 0.05 0.03 0.04 0.06 0.09
2021 0.05 0.03 0.04 0.06 0.08
2022 0.05 0.02 0.04 0.06 0.08
Lasiurus cinereus 2016 0.82 0.71 0.79 0.86 0.91
2017 0.69 0.54 0.64 0.74 0.82
2018 0.73 0.63 0.70 0.76 0.82
2019 0.68 0.59 0.65 0.71 0.77
2020 0.88 0.82 0.86 0.89 0.92
2021 0.79 0.73 0.77 0.81 0.85
2022 0.81 0.75 0.79 0.83 0.87
Lasionycteris noctivagans 2016 0.90 0.82 0.88 0.93 0.96
2017 0.96 0.90 0.95 0.98 1.00
2018 0.97 0.92 0.95 0.98 1.00
2019 0.96 0.91 0.94 0.97 0.99
2020 0.94 0.89 0.92 0.95 0.97
2021 0.86 0.81 0.85 0.88 0.91
2022 0.95 0.91 0.93 0.96 0.98
Myotis californicus 2016 0.96 0.90 0.95 0.97 0.99
2017 0.91 0.81 0.88 0.94 0.97
2018 0.89 0.79 0.86 0.92 0.96
2019 0.95 0.89 0.94 0.97 0.99
2020 0.89 0.81 0.87 0.92 0.95
2021 0.90 0.82 0.88 0.93 0.96
2022 0.79 0.67 0.75 0.82 0.88
Myotis ciliolabrum 2016 0.52 0.33 0.45 0.59 0.72
2017 0.42 0.24 0.35 0.49 0.62
2018 0.57 0.42 0.52 0.63 0.72
2019 0.41 0.29 0.37 0.46 0.55
2020 0.31 0.21 0.27 0.35 0.42
2021 0.41 0.30 0.37 0.44 0.52
2022 0.30 0.20 0.26 0.33 0.41
Myotis evotis 2016 0.88 0.78 0.85 0.91 0.95
2017 0.84 0.71 0.80 0.88 0.93
2018 0.83 0.75 0.81 0.86 0.91
2019 0.84 0.77 0.82 0.86 0.90
2020 0.82 0.75 0.80 0.84 0.88
2021 0.79 0.71 0.76 0.81 0.85
2022 0.76 0.68 0.74 0.79 0.84
Myotis lucifugus 2016 0.84 0.70 0.79 0.89 0.95
2017 0.92 0.82 0.89 0.96 0.99
2018 0.92 0.82 0.88 0.95 0.98
2019 0.94 0.87 0.92 0.97 0.99
2020 0.94 0.86 0.91 0.97 0.99
2021 0.97 0.92 0.96 0.99 1.00
2022 0.97 0.93 0.96 0.99 1.00
Myotis thysanodes 2016 0.57 0.37 0.49 0.64 0.79
2017 0.48 0.30 0.42 0.54 0.66
2018 0.50 0.37 0.45 0.54 0.64
2019 0.32 0.22 0.28 0.35 0.43
2020 0.30 0.22 0.27 0.33 0.39
2021 0.26 0.19 0.23 0.29 0.35
2022 0.27 0.19 0.24 0.29 0.35
Myotis volans 2016 0.78 0.60 0.72 0.85 0.95
2017 0.83 0.65 0.77 0.89 0.98
2018 0.69 0.52 0.62 0.75 0.87
2019 0.53 0.38 0.47 0.58 0.69
2020 0.58 0.46 0.54 0.62 0.71
2021 0.55 0.43 0.51 0.59 0.68
2022 0.54 0.43 0.50 0.58 0.67
Myotis yumanensis 2016 0.63 0.51 0.59 0.67 0.75
2017 0.62 0.51 0.59 0.66 0.74
2018 0.64 0.55 0.61 0.67 0.73
2019 0.59 0.50 0.56 0.62 0.68
2020 0.62 0.54 0.59 0.65 0.70
2021 0.61 0.53 0.58 0.64 0.69
2022 0.57 0.50 0.55 0.60 0.66
Parastrellus hesperus 2016 0.07 0.03 0.05 0.09 0.15
2017 0.05 0.01 0.03 0.07 0.12
2018 0.03 0.01 0.02 0.04 0.07
2019 0.04 0.01 0.03 0.05 0.08
2020 0.04 0.01 0.03 0.05 0.08
2021 0.04 0.01 0.02 0.05 0.07
2022 0.03 0.01 0.02 0.04 0.06
Code
psi_summ %>% 
  filter(type == "Hat", analysis == "OR|WA Only") %>% 
  select(SPECIES, year, mean, q2.5, q25, q75, q97.5) %>% 
  mutate(year = as.factor(year)) %>% 
  gt(groupname_col = "SPECIES") %>%
  cols_label("mean" ~ "Mean",
             "year" ~ "Year",
             "q2.5" ~ "2.5%",
             "q25" ~ "25%",
             "q75" ~ "75%",
             "q97.5" ~ "97.5%") %>% 
  tab_spanner(label = "Quantile", columns = 4:7) %>% 
  fmt_number(decimals = 2) %>% 
  tab_options(row_group.as_column = T)
Table 5: Region-wide summary of the posterior distribution of \(\hat{\psi_t}\) with credible intervals for the analysis including Oregon and Washington only
Year Mean Quantile
2.5% 25% 75% 97.5%
Antrozous pallidus 2016 0.18 0.08 0.13 0.22 0.34
2017 0.18 0.08 0.13 0.21 0.32
2018 0.12 0.05 0.09 0.15 0.21
2019 0.11 0.05 0.08 0.13 0.19
2020 0.09 0.04 0.07 0.11 0.17
2021 0.09 0.04 0.07 0.11 0.16
2022 0.11 0.06 0.09 0.13 0.19
Corynorhinus townsendii 2016 0.12 0.03 0.07 0.15 0.30
2017 0.20 0.06 0.13 0.24 0.43
2018 0.19 0.09 0.15 0.24 0.35
2019 0.20 0.10 0.16 0.24 0.36
2020 0.17 0.08 0.13 0.20 0.30
2021 0.14 0.07 0.11 0.17 0.25
2022 0.13 0.06 0.10 0.15 0.23
Eptesicus fuscus 2016 0.80 0.68 0.77 0.84 0.90
2017 0.82 0.71 0.78 0.85 0.91
2018 0.81 0.72 0.78 0.83 0.88
2019 0.67 0.57 0.63 0.70 0.77
2020 0.71 0.63 0.68 0.74 0.79
2021 0.73 0.64 0.70 0.76 0.82
2022 0.65 0.56 0.61 0.68 0.74
Euderma maculatum 2016 0.02 0.00 0.01 0.03 0.06
2017 0.07 0.03 0.05 0.09 0.13
2018 0.06 0.03 0.05 0.07 0.11
2019 0.06 0.03 0.04 0.07 0.10
2020 0.05 0.02 0.03 0.06 0.08
2021 0.05 0.02 0.04 0.06 0.09
2022 0.05 0.02 0.03 0.06 0.08
Lasiurus cinereus 2016 0.83 0.72 0.80 0.87 0.92
2017 0.73 0.58 0.68 0.78 0.87
2018 0.75 0.65 0.72 0.78 0.84
2019 0.69 0.60 0.66 0.73 0.78
2020 0.86 0.79 0.84 0.88 0.92
2021 0.77 0.69 0.74 0.80 0.84
2022 0.81 0.73 0.78 0.83 0.87
Lasionycteris noctivagans 2016 0.90 0.80 0.87 0.93 0.97
2017 0.96 0.89 0.95 0.98 1.00
2018 0.97 0.91 0.95 0.98 1.00
2019 0.96 0.92 0.95 0.97 0.99
2020 0.96 0.91 0.95 0.97 0.99
2021 0.83 0.76 0.80 0.85 0.89
2022 0.96 0.91 0.95 0.97 0.99
Myotis californicus 2016 0.97 0.93 0.97 0.99 1.00
2017 0.95 0.86 0.93 0.97 0.99
2018 0.94 0.85 0.92 0.97 0.99
2019 0.98 0.94 0.97 0.99 1.00
2020 0.98 0.93 0.97 0.99 1.00
2021 0.96 0.89 0.94 0.98 0.99
2022 0.93 0.84 0.91 0.96 0.98
Myotis ciliolabrum 2016 0.39 0.21 0.32 0.45 0.60
2017 0.28 0.12 0.21 0.34 0.51
2018 0.41 0.25 0.35 0.47 0.58
2019 0.27 0.16 0.23 0.31 0.39
2020 0.23 0.12 0.18 0.27 0.36
2021 0.25 0.15 0.21 0.29 0.37
2022 0.18 0.09 0.14 0.21 0.29
Myotis evotis 2016 0.86 0.75 0.82 0.89 0.94
2017 0.81 0.66 0.77 0.86 0.92
2018 0.81 0.72 0.78 0.84 0.89
2019 0.81 0.73 0.79 0.84 0.89
2020 0.79 0.71 0.77 0.82 0.87
2021 0.74 0.65 0.71 0.77 0.83
2022 0.74 0.65 0.71 0.77 0.83
Myotis lucifugus 2016 0.85 0.73 0.82 0.89 0.95
2017 0.93 0.84 0.91 0.96 0.99
2018 0.92 0.85 0.90 0.94 0.98
2019 0.94 0.88 0.93 0.96 0.99
2020 0.94 0.88 0.92 0.96 0.99
2021 0.96 0.91 0.95 0.98 0.99
2022 0.96 0.91 0.95 0.98 0.99
Myotis thysanodes 2016 0.57 0.36 0.49 0.64 0.80
2017 0.47 0.30 0.41 0.53 0.67
2018 0.49 0.35 0.44 0.54 0.65
2019 0.30 0.20 0.26 0.34 0.42
2020 0.28 0.19 0.24 0.32 0.39
2021 0.23 0.14 0.19 0.26 0.33
2022 0.24 0.15 0.20 0.26 0.33
Myotis volans 2016 0.75 0.57 0.69 0.81 0.92
2017 0.80 0.64 0.75 0.86 0.95
2018 0.64 0.48 0.59 0.70 0.81
2019 0.49 0.36 0.44 0.54 0.65
2020 0.58 0.46 0.54 0.63 0.71
2021 0.51 0.37 0.46 0.56 0.65
2022 0.51 0.38 0.47 0.56 0.65
Myotis yumanensis 2016 0.71 0.59 0.67 0.75 0.81
2017 0.71 0.60 0.67 0.75 0.81
2018 0.73 0.64 0.70 0.76 0.82
2019 0.69 0.60 0.66 0.73 0.79
2020 0.72 0.62 0.69 0.75 0.82
2021 0.68 0.57 0.64 0.72 0.78
2022 0.63 0.54 0.60 0.66 0.73
Parastrellus hesperus 2016 0.06 0.02 0.04 0.08 0.14
2017 0.04 0.01 0.02 0.05 0.11
2018 0.02 0.00 0.01 0.03 0.06
2019 0.03 0.01 0.02 0.04 0.07
2020 0.03 0.01 0.02 0.04 0.07
2021 0.02 0.00 0.01 0.03 0.06
2022 0.02 0.00 0.01 0.03 0.06

JAGS Code

Code

model {
  # Specify priors
  ## Dynamics Parameters
  for (t in 1:(n_years-1)){
    phi[t] ~ dnorm(0, 0.1) #survival
    gamma[t] ~ dnorm(0, 0.1) #colonization
  }

  ## Occurrence-level Parameters
  alpha01 ~ dnorm(0, 0.1) #Year 1 occurrence model intercept

  for(n in 1:n_xcovs){
    alphas[n] ~ dnorm(0, 0.1) #Occurrence covariate intercepts
  }

  ## Detection-level Parameters
  beta0 ~ dnorm(0,0.1) #beta0
  beta1 ~ dnorm(0,0.1) #tmin covariate
  beta2 ~ dnorm(0,0.1) #daylight covariate
  beta3 ~ dnorm(0,0.1) #clut0 covariate
  beta4 ~ dnorm(0,0.1) #Clut1 covariate
  beta5 ~ dnorm(0,0.1) #Clut2 covariate
  beta6 ~ dnorm(0,0.1) #Clut3 covariate (ref is 4)
  beta7 ~ dnorm(0,0.1) #water indicator covariate

  # Fit our Occurrence-level submodel

  ## First Year
  logit_psi[1:n_sites, 1] <- alpha01 + xmat %*% alphas

  ## Following years
  for (t in 1 : (n_years - 1)) {
      logit_psi[1:n_sites, t + 1] <- gamma[t] + xmat %*% alphas #Get gamma colonization term as intercept in following years
  }

  # Ecological submodel: Define state conditional on parameters
  for (i in 1:n_sites){

    z[i,1] ~ dbern(psi[i,1]) #latent occupancy state in year 1 ~ bernoulli with prob based on our covariates and alpha01
    psi[i,1] <- ilogit(logit_psi[i,1])

    #Detection Model Year 1
    for(j in 1:n_visits[i, 1]){
      dets[i, j, 1] ~ dbern(z[i, 1] * p[i, j, 1])
      p[i, j, 1] <- ilogit(beta0 +
        beta1 * tmin[i, j, 1] +
        beta2 * dayl[i, j, 1] +
        beta3 * clut0[i, j, 1] +
        beta4 * clut1[i, j, 1] +
        beta5 * clut2[i, j, 1] +
        beta6 * clut3[i, j, 1] +
        beta7 * wind[i, j, 1])
    }
    
    #Occurrence Probs and Latent state for years 2 to n_years
    for (t in 2:n_years){
      z[i,t] ~ dbern(psi[i,t])
      psi[i,t] <- ilogit(logit_psi[i, t] + z[i, (t-1)] * phi[t-1]) 
    
      #Detection model for years 2 to n_years
      for(j in 1:n_visits[i, t]){
        dets[i, j, t] ~ dbern(z[i, t] * p[i, j, t])
        p[i, j, t] <- ilogit(beta0 + 
        beta1 * tmin[i, j, t] +
        beta2 * dayl[i, j, t] +
        beta3 * clut0[i, j, t] +
        beta4 * clut1[i, j, t] +
        beta5 * clut2[i, j, t] +
        beta6 * clut3[i, j, t] +
        beta7 * wind[i, j, t])
      }
    }
  }
  # Derived parameters
  # lambda measurements (trend)
  lam.tot<-psi.hat[7]/psi.hat[1] 
  lam.avg<-mean(psi.hat[])
  lam.tot.avg<-avg.psi[7]/avg.psi[1]
  lam.avg.avg<-mean(avg.psi[])
  
  #Original estimate for psi.hat (just the intercept of our original model)
  psi.hat[1]<-ilogit(alpha01)
  
  #Psi estimates for following years
  for (t in 2:n_years){
    col[t]<-ilogit(gamma[t-1])
    surv[t]<-ilogit(gamma[t-1]+phi[t-1])
    psi.hat[t]<-psi.hat[t-1]*surv[t]+(1-psi.hat[t-1])*col[t]
    turnover[t]<-((1-psi.hat[t-1])*col[t])/((1-psi.hat[t-1])*col[t]+surv[t]*psi.hat[t-1])
    }
  #Avg psi at year 1
  avg.psi[1] <- mean(psi[ ,1])
  
  #Avg psi at year t
  for(t in 2:n_years){
    avg.psi[t] <- mean(psi[,t])
  }
    
  for (t in 1:n_years){
    p.hat[t]<-ilogit(beta0)
    }
}

References

Kéry, Marc, and J. Andrew Royle. 2021. “Chapter 4 - Modeling Species Distribution and Range Dynamics, and Population Dynamics Using Dynamic Occupancy Models.” In, edited by Marc Kéry and J. Andrew Royle, 207–97. Academic Press. https://doi.org/https://doi.org/10.1016/B978-0-12-809585-0.00004-1.
Reichert, Brian E., Mylea Bayless, Tina L. Cheng, Jeremy T. H. Coleman, Charles M. Francis, Winifred F. Frick, Benjamin S. Gotthold, et al. 2021. “NABat: A Top-down, Bottom-up Solution to Collaborative Continental-Scale Monitoring.” Ambio 50 (4): 901–13. https://doi.org/10.1007/s13280-020-01411-y.
Rodhouse, Thomas J., Rogelio M. Rodriguez, Katharine M. Banner, Patricia C. Ormsbee, Jenny Barnett, and Kathryn M. Irvine. 2019. “Evidence of Region-Wide Bat Population Decline from Long-Term Monitoring and Bayesian Occupancy Models with Empirically Informed Priors.” Ecology and Evolution 9 (19): 11078–88. https://doi.org/10.1002/ece3.5612.
Rodhouse, Thomas, and Kathryn Irvine. 2017. The Master Sample Concept as a Catalyst for Collaborative Continental-Scale Research and Monitoring.
Udell, Bradley J, Bethany R Straw, Tina Cheng, Kyle D Enns, Winfred Frick, Benjamin Gotthold, Kathryn M Irvine, et al. n.d. “Status and Trends of North American Bats Summer Occupancy Analysis 2010-2019 Data Release.” https://doi.org/10.5066/P92JGACB.